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ABSTRACT 


This thesis will examine the effect of nonlinearities on the propagation of orbit 
uncertainties in order to gain insight into the accurateness of the estimation of covariance 
with time. Many real-world applications rely on a first-order approximation of nonlinear 
equations of motion for propagation of orbit uncertainty. The nonlinear effects that are 
ignored during the linearization process can greatly influence the accuracy of the 
solution. A comparative analysis of linear and nonlinear orbit uncertainty propagation is 
presented in order to attempt to determine when linearized uncertainty becomes non- 
Gaussian. An examination of performance metrics is then accomplished to compare 
linearly propagated uncertainty to uncertainty propagated using a second-order 
approximation. An attempt is then made to develop a performance metric that determines 
when propagated uncertainty is no longer Gaussian. The results show it is difficult to 
determine a clear method of defining when the linear approximated uncertainty is no 
longer Gaussian, but there are metrics that can be implemented given a user-defined 
threshold of performance. 
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I. INTRODUCTION 


This thesis will examine the effeet of nonlinearities on the propagation of orbit 
uneertainties in order to gain insight into the aeeurateness of eovarianee estimation with 
time. Aeeurate eovarianee ean be a valuable tool in spaee situational awareness. 
Covarianee information ean be used in the eorrelation of uneorrelated traeks (UCTs) 
[l]-[3], eomputing probability of eollision [4] and sensor tasking [5]. 

Air Foree Spaee Command (AFSPC) eurrently traeks over 22,000 orbiting objeets 
via the Air Foree Spaee Surveillanee Network (SSN). When the SSN deteets an objeet 
that does not eorrelate to an objeet in the Spaee Objeet Catalog, it is eonsidered a UCT. 
Dr. Alfriend proves in [1] that the use of eovarianee for uneorrelated traek assoeiation is 
the optimal measure as long as the uneertainties remain Gaussian. Provided the 
eovarianee is aeeurate, it ean then be used to assoeiate the UCT in a more statistieally 
valid and automated manner than the eurrent labor-intensive proeess [2]. In addition to 
UCT assoeiation, eovarianee is also useful in ealeulating the probability of eollision 
between orbiting objeets. 

The inereased relianee on spaee systems in reeent history has led to inereasing 
probability of eollisions of on-orbit assets. A less aeeurate method of the eollision 
avoidanee ealeulation involves predieting eonjunetions that fall within boxes eentered on 
the two orbiting objeets in question. The disadvantage of this is that it does not take into 
eonsideration the aeeuraey of the predieted boxes. Prior to launeh of the International 
Spaee Station, NASA began using eovarianee for the probability of eollision ealeulations. 
Other Department of Defense and national ageneies have begun to use the same approaeh 
for predieting eonjunetions and the ensuing probability of eollision of on-orbit assets. 
Mueh like UCT assoeiation, an aeeurate eovarianee is a valuable tool in determining with 
greater preeision the probability of eollision between two objeets that will pass elose to 
eaeh other [4]. In terms of determining whether spaeeeraft eollision avoidanee 
maneuvers are neeessary, the aeeuraey of the probability of eollision ealeulation ean then 
be a key faetor with respeet to propellant eonservation and on-orbit safety. 
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Most applications involving uncertainty propagation involve a linear 
approximation of a mathematical model of the nonlinear equations of motion. The 
linearization of the equations of motion for uneertainty propagation simplifies things a 
great deal; however the nonlinear effeets that are ignored during this linearization proeess 
ean greatly influence the aecurateness of the solution. With this in mind, this thesis 
primarily attempts to answer two main questions. First, ean examining the performance 
of nonlinear uneertainty propagation provide any insight into when the linear 
approximation fails? Second, which errors in the initial orbit uncertainty cause the 
propagated covariance to be non-Gaussian and how large do these errors have to be? 

Through comparative analysis of linear and nonlinear orbit uncertainty 
propagation techniques, we ean attempt to determine when the effects of the 
nonlinearities eause the inaccuracy of the covariance to be large enough to render it no 
longer useful. In order to gain insight into when the linear approximation fails, this thesis 
will consider the elassical orbital elements and two-body motion with a comparison of a 
first-order approximation to a second-order approximation. 

Chapter II is a baekground ehapter. It gives a brief review of the Clohessy- 
Wiltshire equations, the definition of the orbit state that will be used and then an 
overview of the higher order approximation that will be used in the comparative analysis. 
Chapter III presents the results of the eomparative analysis and Chapter IV addresses the 
topic of determining appropriate metrics to assess the performance of a first-order 
approximation. 
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II. BACKGROUND 


A, FORMULATION OF THE PROBLEM 

1, Overview 

Expressions that describe the dynamics of orbital motion are based on a set of 
nonlinear differential equations. However, in many real-world applications, the 
propagation of uncertainty is accomplished using a first-order approximation. Given 
initial Gaussian uncertainty, the propagated uncertainty remains Gaussian as long as the 
system remains linear. The system however is nonlinear, and eventually the neglected 
nonlinearities influence the accuracy of the propagated covariance. 

The dynamics of uncertainty propagation are often described using a first-order 
approximation. More complex models exist; however selecting the appropriate model is 
a balance between simplicity and accuracy. With increased accuracy, comes higher cost 
in terms of computation time, computing resources and overall problem complexity. 
Insight into when a particular mathematical model is no longer valid or when it no longer 
meets requirements is very useful. Determining this threshold will be the primary focus 
of this thesis. Previous work by Fujimoto et al. [6], and Park and Scheeres [7], has 
shown significant increases in model accuracy can be obtained by using a second-order 
approximation. Much of the comparative analysis presented in this thesis relies on a 
comparison of second-order to first-order approximation for covariance propagation. 

The background information presented in the remainder of this chapter consists of 
a brief review of the Clohessy-Wiltshire (CW) equations, the definition of the orbit state, 
and then presents an overview of the higher order approximation that will be used for 
comparative analysis. 

2. Linearized Model 

A linearized mathematical model is useful for describing the motion of an orbiting 
vehicle relative to a reference orbit. The well-known Clohessy-Wiltshire equations 
provide a linearized model that approximates relative orbital motion. The CW equations 
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are the linearization of the equations of motion about a eireular orbit and spherieal Earth. 
Although there are more eomplieated linearized models that eonsider non-eireular orbits 
and perturbations [8], the CW equations ean be useful for demonstrating the types of 
relative motion obtained by elosely orbiting objeets. The derivations of the CW 
equations ean be found in found in [9]. For purposes of this thesis, the CW equations are 
merely stated. The CW equations are: 

x-2ny-3n^x = 0 

y + 2nx = 0 (1-1) 

z + n^z = 0 

where n is the mean motion defined by: 


n = 



( 1 . 2 ) 


The solutions for the CW equations are given by the following: 




yn j 


X = (4-3oos(nt))xQ + — |sm(nt) + j |(^l-oos(nt)) 


y = 6 (sin (nt) - nt^ Xq+2 


(x \ 
-^0 

\n ) 


cos 


(nt) + yo+ ^j(4sm(nt)-3«t) (1.3) 


z = Zq cos(nt) + | — |sm(nt) 
n 


The CW frame is illustrated in Figure 1. The x-axis or radial component lies 
along the object’s radius vector. The y-axis or in-track component lies in the orbit plane 
perpendicular to the x-axis. The z-axis or cross-track component is perpendicular to the 
orbit plane. The origin of the CW frame moves with the reference object and is a rotating 
frame. The angular velocity of the rotating frame is equal to the orbital angular velocity. 
This coordinate system is also referred to as the radial, in-track, cross-track or RIC frame 
and is very useful when discussing differences between two orbits. 
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The motion of an anecdotal low earth orbit object under the CW model is 
illustrated in Figure 2. One can see that the motion in the radial and cross-track 
directions are periodic and do not grow in magnitude under the linear assumption. The 
in-track motion does change significantly with time. It is the growth in the in-track 
direction that has higher order effects on the linear model and causes the uncertainty to 
no longer be Gaussian. The goal is to determine when this happens as a function of time 
and initial uncertainties. The secular term from Equation (1.3) is given by: 

Secular Motion = (-bnXj, (1-4) 



In terms of differential motion in the orbital elements in the CW model, 
inclination, right ascension of the ascending node and eccentricity all cause periodic 
variations. The in-track secular motion is caused by differentials in the mean anomaly 
and semi-major axis. In an attempt to isolate the non-Gaussian behavior caused by the 
secular motion in the in-track direction, the two-dimensional state of mean anomaly and 
semi-major axis will be used. This 2D state should capture the essential elements of the 
nonlinear behavior that causes the non-Gaussian behavior of the uncertainty. 
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Motion described by Clohessy-Wiitshire Equations 


X 10* X 10* 




Orbits 



Figure 2. Motion Described by CW Equations 


3, Definition of the State 

In an attempt to isolate the non-Gaussian behavior caused by the secular motion 
in the in-track direction, the state considered for this thesis will consist of the semi-major 
axis and mean anomaly. The dynamics of an n dimensional system can be expressed as: 

(1.5) 



where the solution to Equation (1.5) is given by: 

x(t) = (^(t,Xo,to) (1.6) 


and Xq is the initial state at time t^. Given the reference solution x(t,to), the state 
deviation vector can be represented by: 

Sx = x[t)-x[t) (1.7) 


It follows then that: 

Sx = f{t,x)-f{t,x) (1.8) 


The analysis presented will consist of the two dimensional state consisting of the 
change in semi-major axis and mean anomaly. The “truth” state consists of two-body 
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motion with no perturbations. The equations of motion for the two dimensional state are 
expressed in Equation (1.9): 


a = ag 
d = 0 



(1.9) 


The highest order approximation used for this analysis will be to seeond-order. 
Although a higher order approximation does provide a more aeeurate solution, the most 
signifieant performanee gain over the linear approximations eomes from the seeond-order 
approximation. To obtain the equations that approximate the motion, we will expand in a 
Taylor series about the reference orbit to second-order. 


The state is given by: 


Sx = 


5a 

5M 


( 1 . 10 ) 


where the notation indicates a normalized dimensionless term. The non-dimensional 
change in semi-major axis is given by: 



a 


( 1 . 11 ) 


where the “bar” notation denotes the value on the reference orbit. The time scale used 
will be expressed in terms of number of orbits, r. The time scale is given by: 


T 


nt 

Ik 


( 1 . 12 ) 


The expressions for the change in semi-major axis and mean anomaly (expanded to 
second-order) become: 
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^ * o * 

oa -oar. 


5M = 


« e * \S7t / ^ 
-inoa H- yoa J 


t + SM. 


(1.13) 


Finally, using the first-order term from Equation (1.13), the linearized state transition 
matrix is: 


®(r) 


^ 1 0^ 
\j 


(1.14) 


4. Nonlinear Uncertainty Mapping 


Equation (1.14) is the familiar linear state transition matrix determined by a first- 
order approximation. In order to eapture the nonlinear behavior negleeted in the linear 
approximation. Park and Seheeres [10] proposed a higher order method of approximating 
an orbit state and eovarianee by ineorporating the higher order Taylor series terms. The 
higher order state solutions as a function of initial conditions are referred to as state 
transition tensors (STT). 


Erom Park and Seheeres [10], the state deviation vector the mean deviation 
vector 5m, as well as the uncertainty P, expanded to second-order are given by: 


= (1). d). P 

T i,aT j,a a 


0 

aa 


- 5m,5m, + \<t>i,aAaP + P'pPL) 


(1.15) 

(1.16) 
(1.17) 


Where 5xg and P” are the initial state and covariance, respectively. Einstein’s 
summation notation is used here and the subscripts and “k”, respectively represent the 
ith and kth component. Given the state vector from Equation (1.10) and (1.13), the state 
transition tensors expressed non-dimensionally are: 
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® u .(0 = 


dSx 

dSx 

d^Sx 


y-lnr 1 j 




ddxdda 
d^dx 


^ 0 
15 

V 


0 


m 0 
2 y 


ddxddM 


0 0 
vO Oy 


(1.18) 


The first of Equation (1.18) represents the familiar linear state transition matrix. 
The second and third of Equation (1.18) are the STTs and represent the contributions 
from the second-order Taylor series expansion. One can see from the expressions for the 
second-order STTs in Equation (1.18), all the higher order terms are zero with the 
exception of the Ojn term, which represents the second partial derivative of SM with 
respect to ((5'a*] . 


d) 


2,11 


d^SM 

dSa^ 


(1.19) 


Due to many of the higher order terms in the STTs being equal to zero, the 
expressions for the state, mean, and particularly the covariance can be simplified a great 
deal. The expression for the state given in Equation (1.15) expanded from Einstein’s 
summation notation is written as: 


5x = 


^Sx^ + (j)^ ^ 5 x\ + ^ 1 ^^x^^x^^ + (j)^ ^ 2 ^x^Sx 2 + (j)^ 21^^2^22^^2 ) 
^ 2^2 iSx^ +^2 2^^2 l^X^Sx^ + ^ 2^2 i2Sx'^SX2 + ^ 2^2 21^^25X° + ^ 2^2 22 ^ ^2^ A ) 


( 1 . 20 ) 


Given initial Sx^: 


Sx„ = 


Sal 

SM, 


( 1 . 21 ) 
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when one simplifies Equation (1.20) and diseards the terms that equal zero, the equation 
for the state deviation veetor as a function of non-dimensional time r is: 




( 

-l>7rSal + 

V 


Sal 






t + SM, 


( 1 . 22 ) 


Note that the second term of Equation (1.22) is simply the second-order Taylor 
series expansion of SM shown in Equation (1.13). In fact Equation (1.22) and Equation 
(1.13) are equivalent expressions. 

The second-order mean deviation vector from Equation (1.16) is expressed as: 


(5'm. = 



■*" + ^ 1 , 21 -^,! + 22 ^ 2,2 ) 
^ 2 ,\ 2^\,2 ^ 2 , 21 ^ 2 ,\ ^ 2 , 22 ^ 2,2 ) 


(1.23) 


however, once simplified by combining with Equation (1.18), the mean shown in 
Equation (1.23) reduces to: 


(5'm(r) 


15;r 


P^T 


(1.24) 


The expansion from Einstein’s notation of the second-order covariance is left for 
the Appendix. After combining the expression for the covariance in Equation (1.17) and 
the STTs from Equation (1.18), and simplifying terms equal to zero, the approximation 
for the second-order covariance can be expressed as: 


P = 

y 



KA^Pn 


j) 

riy 11 


+^2,2^2, 


2-* 22 


Sm^Sm^ 


1 

+ — 
4 


> 2 , u 4 u (« 


+ 

^ ll-* 11 


+ /? 


11^11 ) _ 


(1.25) 


or as a function of time t , 
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P(r) 


-'ll 





(1.26) 


Using Equation (1.22) and (1.26) we can now propagate the state and covariance 
using a second-order approximation. 
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III. COVARIANCE ANALYSIS 


A, INTRODUCTION 

This section examines covariance propagation through comparison of uncertainty 
ellipsoid propagation and Monte Carlo simulations. The cumulative probability 
distribution for first-order and second-order approximations compared to the theoretical 
distribution for a linear system is also presented. To begin, we will examine how 
uncertainty propagates under the first-order and second-order approximations. 

B, ELLIPSE PROPAGATION 


For the uncertainty ellipse propagation, a demonstration of how an initial 
Gaussian distribution propagates under unperturbed two-body equations of motion is 
presented. Only the two dimensional state consisting oiSa* and SM is considered. 

For the error ellipsoid analysis presented in this section, 10^ points on the initial 
3cr uncertainty ellipse are propagated using the linear approximation given by the state 
transition matrix in Equation (1.14). The points are also propagated using the second- 
order approximation given in Equation (1.13). A Monte Carlo simulation is then 
performed by propagating 10^ normally distributed points obtained from the initial 
uncertainty. The Monte Carlo points are propagated using the complete system 
dynamics: 


™ ™ « 
da = oa^ 


( 


dM=M-M = 


A 


(a +da^ 


+Mq 


+ Mq 

a 




( 2 . 1 ) 


Note that Equation (2.1) is not expressed in dimensionless terms or in the 
timescale r. The following shows the expressions for the state in non-dimensional terms. 
The non-dimensional semi-major axis is: 
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* Cl 

a = — 
a 


where the “bar” over the variable represents a value on the referenee orbit, 
dimensional time is given by: 


( 2 . 2 ) 


The non- 


T = ■ 


nt 

In 


(2.3) 


Given that the semi-major axis is: 


a = a + da 


(2.4) 


the normalized semi-major axis is 


a =\ + Sa 


(2.5) 


The expression for mean anomaly is: 


M = nt + Mr, 


( 2 . 6 ) 


Substituting z into Equation (2.6) for non-dimensional terms yields: 


M = 2n 


a' 

— r + Mq = 2nA^T + 
\n ) \ a 

3/2 


( *\-M 2 

a j T + 


(2.7) 


Note that: 


dM 

dr 


dz 

= 2n(ay'" 


( 2 . 8 ) 


Therefore, given that 5M 


M -M the exaet solution for the ehange in mean anomaly is: 




(a*p-l 


Z + dMr, 


(2.9) 
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1, Ellipse Propagation Results 


By propagating the 3fJ uncertainty ellipsoid using first and second-order 
propagations and comparing to the exact propagation of the Monte Carlo points, we are 
able to examine the validity of the approximations at various propagation times. For this 
scenario, the uncertainty is given by: 



( 2 . 10 ) 


where the semi-major axis uncertainty is 21 km, the mean anomaly uncertainty is 0.01° 
and the reference semi-major axis is a = 7000 km. Given these parameters, the orbital 
period is 1.619 hours. The standard deviations are given by: 



The initial distribution as well as the 3cr ellipse based on the initial uncertainty 
are shown in the top left of Figure 3. The remaining plots are the 3cr ellipse propagated 
using first and second-order approximations, as well as the Monte Carlo points 
propagated using the complete system dynamics at 2, 5, and 10 orbits. 

The green ellipse represents the 3cr ellipse propagated using a first-order 
approximation from Equation (1.14) and the blue dashed ellipse is the second-order 
approximation from Equation (1.13). The Sa* -JMplane is rotated so the principal axes 
of the linearly propagated ellipse are aligned with the horizontal and vertical axis of the 
plot. Ideally, the propagated ellipse should continue to capture a large number of the 
Monte Carlo sample points. One can see that the accuracy of the linearly propagated 
ellipse capturing the sample points significantly degrades with time. 
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X Error Distribution and 3a Ellipse 


3a Co\«riance After 2 Orbits 



Figure 3. Initial Gaussian Distribution and Propagated Gaussian Distribution with 

Corresponding 3a Ellipsoids 


The theoretieal distribution for a two-dimensional system is given by; 

F{k) = \-e^ (2.12) 

where k represents the number of standard deviations. For A=3, the theoretieal 
probability distribution is 0.9889, i.e., the 3a ellipse should eapture 98.89% of the Monte 
Carlo sample points. As seen in Table 1, for the initial uneertainty given in (2.11), by 10 
orbits the linear propagation has failed to represent a significant portion (nearly half) of 
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the system dynamics. The illustrations in Figure 3 show the second-order propagation 
does a much better job of characterizing the full system dynamics. The results of the 
probability distribution captured by the first-order 3cr ellipse are summarized in Table 1. 


Theoretical 

Initial 

2 Orbits 

5 Orbits 

10 Orbits 

0.9889 

0.9992 

0.8383 

0.6512 

0.5045 


Table 1. Linear Probability Distribution through 10 Orbits 


3a Covariance After 0.375 Orbits 



3a Covariance After 0.75 Orbits 




Figure 4. Propagated Gaussian Distributions with Corresponding Scr Ellipsoids 
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Although the linear propagation does not fully represent the full system dynamics, 
for a short propagation time, the first-order approximation may be acceptable. One can 
see that at 2 orbits, or 3.24 hours, the first-order approximation still encompasses 83.8 % 
of the sample points. 

Figure 4 illustrates the results of various propagation times out to 1.5 orbits. 
Again, the propagated covariance is rotated so the primary axes of the linear ellipse 
correspond to the axes of the plot. One can see via the distortion of the second-order 
ellipsoid and the Monte Carlo distribution, there are higher order effects even at 0.375 
orbits, or just over 30 minutes of propagation time. Despite some loss from nonlinear 
effects, depending on the accuracy required, the linear approximation at shorter 
propagation times may be sufficient. Table 2 summarizes the numerical results for the 
probability distribution captured by the linear approximation. 


Theoretical 

Initial 

0.375 Orbits 

0.75 Orbits 

1.125 Orbits 

1.5 Orbits 

0.9889 

0.9992 

0.9783 

0.9498 

0.9167 

0.8825 


Table 2. Linear Probability Distribution through 1.5 Orbits 


2, Ellipse Rotation 

The ellipse rotation consists of rotating the major and minor axis of the first-order 
ellipse to correspond to the primary axes of the plot. This allows a closer examination of 
the higher order effects. The second-order contribution to 5M comes from the following 
term of Equation (1.22): 

(2.13) 

Figure 5 shows the non-rotated covariance ellipsoids propagated to two orbits 
using the same standard deviations given in (2.11). The plot on right of Figure 5 is 
zoomed in to show the second-order ellipse curvature. 
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Non-rotated Ellipsoid After 2 Orbits Non-rotated Ellipsoid After 2 Orbits 



Figure 5. Non-Rotated Covariance Ellipsoids after Two Orbits 


The higher order contribution to dM is small compared to the total value of dM 
but it is large compared to the ellipsoid minor axis. A property of linearly propagated 
covariance is that the area of the linear ellipsoid is constant with time. The covariance 
does not grow in the da direction, so as it grows in the SM direction, the ellipse becomes 
thinner with time along the minor axis. 

To examine the higher order effects it is necessary to determine a coordinate 
system where the first-order covariance is diagonal. This will rotate the axes of the 
linearly propagated covariance to the primary axes of the plot. To do so, we use the 
following; 

Jy = RJx (2.14) 


where dx is the state given Equation (1.22) and dy is the rotated state; 


dy = 


ysM 

ysa_ 


(2.15) 


R is the clockwise rotation matrix; 
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R = 


cos 6' 
-sin 6' 


sin 6' 
cos6' 


The covariance rotation is accomplished by; 

Pv = RPxR" 

where is the linearly propagated covariance 



and the rotated covarianee P^, is given by: 


P> = 


P P 


->•12 


>12 


' yl2 


(2.16) 


(2.17) 


(2.18) 


(2.19) 


which gives: 

+ -/:22)':os(2O) + P„jSm(20) 

P,i,=\(P.n +P,22 )-^(Au-P, 22)^os(20)-P,„sm(20) (2.20) 

Py,2=P2i2<x>s(20)-^{P,t,-P,,i)sm{20) 


To obtain the coordinate system where the covariance is diagonal, we set P^^^ = 9 
and solve for the rotation angle 6 ; 


9 = — aretan 
2 


2P 


xl2 


P - P 

V xl 1 x22 y 


( 2 . 21 ) 


or 
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( 2 . 22 ) 


6 = — arctan 
2 


2 ( 

^Sa^^Su) 




The rotation angle obtained from Equation (2.22) is then used with the rotation matrix 
given by Equation (2.16). 

C. PROBABILITY DISTRIBUTION 

1. Overview 


This seetion examines the eumulative probability distribution using the first and 
second-order approximation as well as the exact solution for covariance propagation. 
The results of the probability distribution are compared to the theoretical cumulative 
distribution function to assess the performance of the approximation. The Gaussian 
multivariate probability density function is; 


/(X) 




(2.23) 


where N is the dimension of the state. The probability density function as a function of 
number of standard deviations, k is: 


m= 


1 



(2.24) 


The Gamma function r(x) is given by: 

The probability distribution function as a function of k is then given by: 


(2.25) 
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F{k) = 


1 


(2.26) 


—1 

2^ r 




0 


For our purposes N=2. The theoretical probability density function and probability 
distribution function for a two dimensional system are; 


f{k) = ke^ 


F{k) = l-e 2 


(2.27) 


2. Probability Distribution Results 


For the probability distribution, five evenly spaced points were investigated over a 
range of uncertainty of the normalized semi-major axis. The range used for this analysis 
is 1x10'"^ to 3x10'^ corresponding to a range of semi-major axis uncertainty of 0.7 km to 
21 km. This is summarized in Table 3. Initially the uncertainty in mean anomaly is left 
at 0.01° through the range of semi-major axis uncertainty. 


Range of Uncertainty in Normalized Semi-Major Axis 


o-,. 


1.000x10" 


1.250x10" 


1.550x10' 


2.275x10' 


3.000x10' 


Table 3. Normalized Semi-Major Axis 


Given Equation (2.27), at k=3 for the theoretical linear system, F(k) is 0.9889. 
That is, a 2D system at 3cr standard deviation should theoretically capture 98.9% of the 
distribution. In Table 4, the data presented shows the value of k when 98.9% of the 
distribution is captured for each of the methods of propagation. If the approximation is 
adequately capturing the dynamics, the k @ 98.9% row should remain close to 3. If it 
does not, this indicates more standard deviations are required to capture the system 
dynamics. The values of F(k) at A=3 for the first and second-order approximation as well 
as the exact propagation are is also shown in Table 4. This allows us to evaluate how 
much the propagation degrades with time. 
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During this analysis, it was discovered that the cr^^. =1.0x10^ / =0.0Tease 

showed exeellent agreement in the first-order approximation. This is due to the higher 
order effeets being a fimetion of 5a and the pereentage ehange between the first and 
seeond-order approximations being smaller with larger cr^^. The results presented in this 
seetion start at cr^^. = 8.25x10 or 5.775 km. The behavior eorresponding to the 
(j^^. = 1.0x10 case will be presented later along with diseussion on the relative behavior 
between cr,. and cr.,,. 

Sa oM 



Probability Distribution - 8.25x10 '' 

Orbits 

2 

4 

6 

8 

10 

20 

40 

60 

80 

100 

Linear 

k @ 98.9% 

3.05 

3.21 

3.51 

3.80 

4.19 

6.51 

11.96 

17.56 

23.30 

29.10 

F(k) @ k=2> 

0.987 

0.983 

0.974 

0.964 

0.951 

0.882 

0.770 

0.684 

0.619 

0.569 

2nd 

Order 

k @ 98.9% 

3.04 

3.16 

3.37 

3.54 

3.73 

4.45 

4.91 

5.00 

5.05 

5.07 

F(k) @ k=3 

0.988 

0.985 

0.979 

0.973 

0.967 

0.949 

0.940 

0.937 

0.936 

0.935 

Exact 

k @ 98.9% 

3.06 

3.19 

3.40 

3.59 

3.79 

4.54 

4.95 

5.06 

5.11 

5.13 

F(k) @ k=3 

0.987 

0.984 

0.977 

0.972 

0.966 

0.947 

0.938 

0.934 

0.933 

0.933 


Table 4. Probability Distribution - = 5.775 km 


One ean see in Table 4 that the first-order approximation begins to signifieantly 
degrade by 20 orbits. This is further illustrated in Figure 6. The dashed red line 
represents the theoretieal probability distribution from Equation (2.27). The remaining 
eurves are the probability distribution plotted at 20 orbit inerements out to 100 orbits. 
The linear propagation eurves flatten out with propagation time. This indieates a larger 
standard deviation is needed to eapture the distribution. 
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Linear Propagation 



k 

Exact Propagation 


^ CT,,. = 0.000825 

n sa 



2nd Order Propagation 


CT . = 0.000825 



k 


0 Orbits j 

20 Orbits 
40 Orbits 
60 Orbits 
80 Orbits 
100 Orbits 
Theoretical 


Figure 6. Probability Distribution Comparison, 100 Orbits - = 5.775 km 


Tables 5 through 7 are the summary of the probability distribution for the 
remainder of the range of uncertainty in semi-major axis. The remaining probability 
distribution comparison plots are left for the appendix. It can be seen in Tables 5 through 
7 that as we increase the initial uncertainty in semi-major axis, leaving the initial 
uncertainty in mean anomaly at 0.01°, the accurateness of the first-order approximation 
decreases. 



Probability Distribution - cr^„. = 1.55x10'^ 

Orbits 

2 

4 

6 

8 

10 

20 

40 

60 

80 

100 

Linear 

k @ 98.9% 

3.68 

5.11 

6.85 

8.78 

10.71 

20.59 

40.94 

61.40 

81.69 

102.06 

F(k) @ k=3 

0.969 

0.924 

0.875 

0.830 

0.792 

0.646 

0.497 

0.417 

0.367 

0.331 

2nd 

Order 

k @ 98.9% 

3.44 

4.09 

4.51 

4.73 

4.85 

5.02 

5.09 

5.10 

5.10 

5.10 

F(k) @ k=3 

0.976 

0.959 

0.948 

0.943 

0.941 

0.936 

0.935 

0.935 

0.935 

0.935 

Exact 

k @ 98.9% 

3.49 

4.18 

4.59 

4.79 

4.93 

5.08 

5.13 

5.15 

5.13 

5.12 

F(k) @ k=3 

0.974 

0.956 

0.946 

0.941 

0.939 

0.933 

0.932 

0.932 

0.933 

0.933 


Table 5. Probability Distribution - <7= 10.85 km 
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Probability Distribution - = 2.275x10'^ 

Orbits 

2 

4 

6 

8 

10 

20 

40 

60 

80 

100 

Linear 

k @ 98.9% 

5.37 

9.38 

13.61 

17.82 

22.15 

44.14 

87.93 

131.9 

175.6 

219.1 

F(k) @ k=2> 

0.916 

0.818 

0.740 

0.682 

0.630 

0.482 

0.355 

0.292 

0.254 

0.226 

2nd 

Order 

k@9S.9% 

4.16 

4.77 

4.92 

4.99 

5.03 

5.10 

5.10 

5.11 

5.10 

5.09 

F(k) @ k=2> 

0.957 

0.942 

0.939 

0.937 

0.936 

0.935 

0.935 

0.935 

0.935 

0.935 

Exact 

k @ 98.9% 

4.26 

4.85 

5.00 

5.06 

5.10 

5.14 

5.13 

5.13 

5.14 

5.13 

F(k) @ k=2> 

0.954 

0.940 

0.936 

0.934 

0.933 

0.932 

0.932 

0.933 

0.933 

0.933 


Table 6. Probability Distribution - - 15.925 km 


Table 7 shows that with an initial uncertainty in the semi-major axis equal to 
3x10' after 2 orbits F(k) at A=3 has already degraded to 84% indicating that the linear 
approximation is no longer representing a significant portion of the system dynamics. 
Figure 7 also illustrates this. In addition to the probability distribution plot, Figure 7 
shows the uncertainty ellipsoid propagation for comparison. 



Probability Distribution - cr^^. = 3x10'^ 

Orbits 

2 

4 

6 

8 

10 

20 

40 

60 

80 

100 

Linear 

k @ 98.9% 

8.30 

15.57 

23.07 

30.61 

38.43 

76.44 

152.6 

228.4 

304.2 

379.9 

F(k) @ k=3 

0.842 

0.709 

0.621 

0.556 

0.507 

0.380 

0.271 

0.222 

0.191 

0.173 

2nd 

Order 

k @ 98.9% 

4.70 

4.97 

5.03 

5.07 

5.10 

5.10 

5.09 

5.08 

5.08 

5.08 

F(k) @ k=3 

0.944 

0.938 

0.936 

0.935 

0.935 

0.935 

0.935 

0.935 

0.935 

0.935 

Exact 

k @ 98.9% 

4.74 

5.03 

5.10 

5.14 

5.14 

5.14 

5.13 

5.13 

5.13 

5.13 

F(k) @ k=3 

0.941 

0.935 

0.933 

0.933 

0.932 

0.933 

0.933 

0.933 

0.933 

0.933 


Table 7. Probability Distribution - =21 km 


It is noteworthy to see that throughout the entire range of uncertainty in semi¬ 
major axis, the second-order approximation is consistent out to 100 orbits and closely 
matches the exact probability distribution. This indicates that for the two dimensional 
semi-major axis/mean anomaly state, the second-order approximation is very close to the 
true system dynamics. 
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Linear Propagation 



2nd Order Propagation 

o, =0.003 
5a* 



k 


Exact Propagation 



k 


0 Orbits 
2 Orbits 
4 Orbits 
6 Orbits 
8 Orbits 
10 Orbits 
Ttieoretical 



Figure 7. Probability Distribution Comparison, 10 Orbits and Gaussian Distribution 

Propagated for 2 Orbits - = 21 km 


3, Summary of Probability Distribution 

The probability distribution is an exeellent way to examine how well a model is 
fully capturing the dynamics of a system. The results indicate that with larger initial 
uncertainty in semi-major axis, the first-order approximation fails much quicker than with 
smaller initial uncertainty. 

The results for the =1.0x10 were not presented in this section and are 
included later in a section that evaluates the relative behavior of cr^^. and 
Additionally, only variations in cr^^. were presented. The results discussing the relative 
behavior will show that with a larger initial relative to cr^^., there is less of a 
percentage difference between the first-order and higher order solutions. 
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IV. PERFORMANCE METRICS 


A. INTRODUCTION 

There has been signifieant work done in the area of nonlinear performanee 
metries. This ehapter is an examination of previous work in the area of performanee 
metries of predietion error and uneertainty effeets as well as an examination of the 
relative behavior between initial uneertainty in mean anomaly and semi-major axis. The 
ultimate goal is to determine if a metrie ean be developed to prediet when ignored 
nonlinear effeets on eovarianee propagation have rendered an approximation no longer 
useful. 

B. PREDICTION ERROR COST FUNCTION 
1. Overview 

Horwood, et al. [11] developed a metrie for uneertainty eonsisteney ealled the 
predietion error to help identify errors in nonlinear and linear propagation teehniques. 
The predietion error is given by: 

(x,t )/?2 (x,t)r/x (3.1) 

where and P 2 are the probability density funetions for two satellite states evaluated at a 
common time integrated over the entire state space. The cost metric is defined by: 

c(t) =-In'P£(t) (3.2) 


or 

^(0 = ^(*1 - V (Pi + P 2 )"' (xi - X 2 ) + ^In[det (2;r(P, + P 2 ))] (3.3) 

The first term in Equation (3.3) is one half of the square of the Mahalanobis distance [12] 
and the second term is the determinant of the covariance. It will be useful to analyze 
each of these terms separately: 
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(3.4) 


^l(^)=|(Xl-X2r(Pl+P2r‘(Xi-X2) 
C2(^) = |ln[det(2;r(Pi+P2)) 


where the total cost, c(t) = Cj(t) + C 2 (t). Additionally, for purposes of this thesis, 
although unrealistic, we will be using a “truth” orbit for orbit 2. Therefore, the state 
deviation X 2 and uncertainty P 2 are both zero. Equation (3.4) reduces to; 


.(t) = ^ln[det(2;rPi)] 


(3.5) 


with the sum of Cj and being the entire cost function. The first term of the cost function 
uses the Mahalanobis distance, so prior to evaluating the cost function, we will examine 
the MD with respect to the state. 

2. Mahalanobis Distance 

The Mahalanobis distance is a way to determine the statistically weighted 
difference between two points. The Mahalanobis distance is given by: 

e =5x^V-'dx, 

where (3.6) 

= x(t)-x(t) 


where the “bar” over the state indicates a value on the reference orbit. Recall from 
Equation (1.13) that Sa = Sal reference orbit is the truth orbit. Therefore; 


a-a 


Suq 



SM 


(3.7) 
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Figure 8. Mahalanobis Distance - First and Second-Order Approximations 


The Mahalanobis distance for the first and second-order propagation is plotted in 
Figure 8. For a linear system, the Mahalanobis distance is constant. Our system is 
nonlinear but uses approximations to propagate the uncertainty. As one would expect for 
the initial conditions Sal - ~ ^ the MD starts at 3 and grows larger with 

time. For the linear propagation, the MD starts at 3 and grows infinitely large with time. 
Interestingly, for the second-order approximation, the MD approaches an asymptote. 
How quickly the second-order MD approaches the asymptotic value is proportional to 
Smaller <J^^. results in smaller change in MD suggesting the first and second-order 
approximations show good agreement with the exact propagation, ft will be shown later 
that this is also a function of the relative behavior between initial uncertainty in Sal 
SMq. Overall, the departure of the Mahalanobis distance from the initial value at initial 
time is a method of indicating how well your approximation is capturing the full 
nonlinear dynamics. 
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3. 


Prediction Error Cost Function Results 


The cost function will be used to compare the results of the first and second-order 
approximations. In order to obtain the initial covariance, we must define our initial 
uncertainty. The uncertainty in mean anomaly is left at 0.01°. Equation (1.26) shows 
that it is the uncertainty in the semi-major axis or the radial direction that causes the mean 
anomaly error or in-track error to grow with time. Due to this, we will consider a range 
of uncertainties in semi-major axis. To get an accurate comparison with previous results, 
the non-dimensional semi-major axis uncertainties considered are 1x10'^ to 3x10'^ 
corresponding to a range of uncertainty from 0.7 km to 21 km. One can see from 
Equation (3.3) that the initial cost at time zero is dependent on the selection of the initial 
conditions. The initial condition used for prediction error cost analysis is Sal - 
and SMq =0. A different initial condition would result in a different value for initial and 
final cost and ultimately a different total cost. The choice of initial condition is based 
Sal - being the “worst case” from a range of uncertainty bounded by the 3f7 

ellipse. The total cost for 10 orbits is plotted in Eigure 9. 

Total 1st Older Cost Function fcr 10 Orbits 
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Figure 9. Total Cost Metric for 10 Orbits - First and Second-Order Approximation 
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The results of the total eost as well as the Mahalanobis distanee for the first-order 
approximations are summarized in Table 8. Propagation times are set at 10 as well as 
100 orbits for comparison. Interestingly, for a small initial uncertainty, or cr^^. =1x10 
there is very little change in total cost or MD over both the 10 and 100 orbit propagation 
times. This initially would indicate that the uncertainty consistency is well represented 
for smaller cr^^.. Like the Mahalanobis distance, the performance of the cr^^. = 1x10 case 
is in reality maximizing performance of the relative behavior between and cr^^.. 

The C 2 term of the first-order approximation requires some discussion. The second 
term of the cost is a function of the determinant of the covariance. Alfriend [1] proved 
that for Hamiltonian systems, the volume of the uncertainty ellipsoid is constant with 
time. In other words, the determinant of the covariance is constant. This was also 
demonstrated during the covariance ellipsoid propagation with the “thinning” of the 
ellipsoid as it grew in the 5M direction. Due to the constant area of the covariance 
ellipsoid propagated under the first-order approximation, the difference in the term or 
the total cost in over the propagation time is zero. 


1st Order Approximation 

10 Orbits 


Total Cost 

Term 1 Cost 

Term 2 Cost 

MD 

1.000x10'^ 

1.8440x10'^ 

1.8440x10'^ 

0 

3.000615 

8.250x10'^ 

8.4991 

8.4991 

0 

5.098846 

1.550x10'^ 

105.3632 

105.3632 

0 

14.823173 

2.275x10'^ 

486.5164 

486.5164 

0 

31.337402 

3.000x10'^ 

1463.7658 

1463.7658 

0 

54.189773 

100 Orbits 


Total Cost 

Term 1 Cost 

Term 2 Cost 

MD 

1.000x10'^ 

0.1844 

0.1844 

0 

3.060849 

8.250x10'^ 

849.9114 

849.9114 

0 

41.337911 

1.550x10-^ 

1.0536x10^ 

1.0536x10^ 

0 

145.195203 

2.275x10'^ 

4.8652x10'' 

4.8652x10" 

0 

311.949164 

3.000x10'^ 

1.4638x10^ 

1.4638x10^ 

0 

541.075005 


Table 8. First-Order Approximation Cost and Mahalanobis Distance 


31 





















































The results of the total cost for the second-order approximation are summarized in 
Table 9. As one would expect, overall the total cost for the first-order approximation is 
greater than the second-order approximation indicating greater uncertainty consistency 
with a second-order approximation. However, this does not hold with the cr^^. = 1x10 
case. The total cost for this case is very close for the first and second-order 
approximation. The case with initial uncertainties of cr^^. = 1x10 and = 0.01° show 
good agreement between the first and second-order approximations. This agrees with the 
results obtained from the probability distribution and is also due to the relative behavior 
in initial uncertainties. 

As stated previously, a property of a linearly propagated covariance is that the 
volume is constant with time. The term is simply the volume of the covariance scaled 
by a constant. The condition of the covariance ellipsoid volume being constant with time 
does not hold for the second-order approximation. One can see there is a contribution for 
the Term 2 cost at cr^^. =1.0x10^. This actually results in a slightly higher cost for the 
second-order approximation for this case. 


2nd Order Approximation 

10 Orbits 


Total Cost 

Term 1 Cost 

Term 2 Cost 

MD 

l.OOOxlO'^ 

1.8894x10'^ 

1.8438x10'^ 

4.5560x10'^ 

3.000615 

8.250x10'^ 

6.1524 

5.9763 

0.1761 

4.577401 

1.550x10'^ 

17.7490 

16.8319 

0.9171 

6.531751 

2.275 xlO'^ 

20.7644 

19.1468 

1.6176 

6.877041 

3.000x10'^ 

21.7236 

19.5661 

2.1575 

6.937740 

100 Orbits 


Total Cost 

Term 1 Cost 

Term 2 Cost 

MD 

1.000x10'^ 

0.1873 

0.1827 

0.0045 

3.060305 

8.250x10'^ 

21.5507 

19.6676 

1.8831 

6.952358 

1.550x10-^ 

23.1276 

19.9940 

3.1336 

6.999144 

2.275x10'^ 

23.8234 

19.9231 

3.9003 

6.989003 

3.000x10'^ 

24.2819 

19.8285 

4.4534 

6.975459 


Table 9. Second-Order Approximation Cost and Mahalanobis Distance 


32 
























































4. Summary of Cost Function and Mahalanobis Distance 


Both the prediction error and the Mahalanobis distance give comparison between 
an initial and propagated state. The prediction error takes into account higher order 
statistics with the term. Neither metric provides definitive information or give a 
consistent means of determining when the linear approximation is no longer valid. The 
prediction error cost function is very useful for determining performance gains when 
using higher order filters, but overall doesn’t answer the problem posed in this thesis. 

The Mahalanobis distance gives a statistical difference between an initial point 
and the same point after propagation. Alfriend [1] has shown that the Mahalanobis 
distance can be used for uncorrelated track association in a more statistically valid 
manner. With these methods however, the selection of a specific MD as a metric is a 
balance between accurate correlations and false correlations. Previous work by Hill et al. 
[3] suggest a MD of less than 10 can be used for UCT association with a reasonable 
amount of confidence. There is no definitive total change in MD that identifies when the 
first-order approximation is no longer valid. 

C. UNCERTAINTY METRICS 

In addition to the work with prediction error cost functions, there have been other 
studies into determining metrics to assess and compare performance of higher order 
propagation techniques. In an attempt to determine if these metrics will reveal any 
information into the validity of a first-order approximation, we will examine these 
metrics further. 

Park and Scheeres developed a metric to assess the performance of the State 
Transition Tensor method of mapping uncertainties. It was found that when applied to 
the state used in this thesis, this metric displayed some singularities close to the initial 
time, r = 0. Due to this, it was required to slightly modify the Park and Scheeres metric. 
The results of the modified metric are presented in this section. To begin, a review of the 
Park and Scheeres metric is conducted. 
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1 . 


Park and Scheeres Nonlinearity Metric 


The Park and Scheeres [7] local nonlinearity index was proposed to assist with 
determining the sufficient order of higher-order solutions. This local nonlinearity index 
is defined as; 




m 



sup 


o 

O 



(7x;( 

o 

o 

X 



(3.8) 


where is the mth order approximation of Equation (1.8), dx\ is the kth 
initial condition for (^x chosen from the boundary of the initial confidence region, i.e. on 
the surface of the kcr ellipse, and is the exact solution of Equation (1.8). Einstein’s 
summation notation is used and n is the dimension of the state while N is the number of 
initial conditions. 
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Figure 10. Park and Scheeres Nonlinearity Index 


Upon examination of this metric, it was discovered that for the state used in this 
thesis, there were singularities due to the denominator of Equation (3.8) occasionally 
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being very elose to zero, partieularly near time t = 0 . The singularities are illustrated in 
Figure 10. The results in Figure 10 were obtained using Equation (3.8) with 1000 initial 
conditions sampled on the 3cr ellipsoid with cr^^. = 0.003 and = 0.0T. 


The steady state value of rj for this particular set of initial conditions and initial 
uncertainty is actually 1.1255x10' for the first-order approximation and 1.1823x10' for 
the second-order approximation. One can see the singularities at the beginning of the 
propagation are distorting the curve. The following will derive the limit of the index in 
order to determine what the steady state behavior is. 

Recall that the expression for the exact solution of 5M is; 


* r / * \ “3/2 

5M = In (a ) -1 z + dM, 


/ « \ —3/2 

The Taylor series expansion of (a j is given by: 


(cT-l + tc,{Say 

./=1 

j\ da 


(3.10) 


therefore the exact solution for SM becomes; 


00 

SM* = In'^cj {^Sa*yT + SMq 


(3.11) 


and the mth order approximation for the change in mean anomaly is; 


m 

SW" = 2nY,Cj [Sa*)'z + SM, 


(3.12) 


The c , to m=3 are; 




(3.13) 



Although the nonlinearity index is a function of the entire state, da is constant. 
Therefore, only the expression for SM need be considered for the index. Substituting 
Equation (3.11) and (3.12) into Equation (3.8), the local nonlinearity index becomes; 
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sup 

k=U...,N 


l/TT 

^c.(daiy (dal)' 

7=1 7=1 


00 

2;rY,Cj(dal)'T + dMl 

7=1 



(3.14) 


Consolidating the sigma notation in the numerator gives; 


77 


m 



sup 

k=l,...,N 


27tT 

00 



y=m+i 


00 

2;r^c 

7=1 

■XSa\)r + SMl 


(3.15) 


With the results in Equation (3.15), there are two observations about the 

performance of this metric applied to our state. One is that when SM^ and 

00 

2nY,Cj(Sal)' T are equal and of opposite sign, rj"' will be infinite. These singularities 

7=1 

often appear near r = 0. Due to this it is necessary to slightly modify the local 

nonlinearity index. The second observation is that eventually the j r term will 

0 

dominate dM,^. Once this occurs, we can express the steady state value reached by 
Equation (3.15) as; 


sup 




j=m+l 


k=l,...,N 


7=1 


sup 

k=l,...,N 


''m+l 




(3.16) 
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2. Proposed Performance Index #1 


In order to mitigate the singularity behavior with the Park and Scheeres local 
nonlinearity index, the problem with the term in the denominator being zero needs to be 
addressed. The first of two proposed performance metrics is given by; 


11 
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sup 


5x”' 

O 

O 

)-5x* 

{t;5xlt^) 


5x* ( 

t;5xl/ 

') 



(3.17) 


It is important to note that(7x must be dimensionless or each term in dn must have the 
same dimension. Again, due to the fact that 5a is constant, we can simplify the 
expression for rj and express it in terms of 5M. Equation (3.17) becomes: 


7 




sup 


5M"‘[f,5xlf)-5M*[f,5xlf) 


[5a*)\[5M*)^ 


1/2 


(3.18) 


Equation (3.18) will start at zero and asymptotically approach the following 

value; 
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sup 

k=\,...,N 



(3.19) 


Consider the k<j . to be the maximum value of 5a*. This occurs when 5Mr. = 0. This 

oa 

leaves: 
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m 




(3.20) 


/ »\- 3/2 

where the Cj come from the Taylor series expansion of (a 1 about the reference value 
a* = I given in Equation (3.10). 
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1st Order Case 
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Figure 11. Proposed Performance Metric #1-10 Orbits 


Figure 10 shows Performance Index #1 for the range of cr^^.. The uncertainty in 
mean anomaly is left at 0.01°. One can see that for the (J^^. =1.0x10^ case, the first-order 
is still showing some discontinuities during initial propagation time and settling on the 
theoretical value. It is unclear what is causing these discontinuities. Based on the results 
of this metric, it provides a measure of increased performance with a higher order 
solution, but does not provide information on when the first-order approximation is no 
longer valid. Due to this, further investigation into the discontinuities is not warranted. 



1.000x10'^ 

8.250x10'^ 

1.550x10'’ 

2.275x10'’ 

3.000x10'’ 

1 St Order Theoretical 

3.7500x10'^ 

3.0937x10'’ 

5.8125x10'’ 

8.5312x10'’ 

1.1250x10'’ 

1 St Order Numerical 

3.7518x10-^ 

3.0936x10'’ 

5.8134x10'’ 

8.5341x10'’ 

1.1255x10'’ 

2nd Order Theoretical 

1.3125x10'’ 

8.9332x10'*’ 

3.1532x10'’ 

6.7930x10'’ 

1.1812x10''' 

2nd Order Numerical 

1.3127x10'’ 

8.9338x10'*’ 

3.1543x10'’ 

6.7972x10'’ 

1.1823x10'" 


Table 10. Proposed Performance Metric #1-10 Orbits 
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3. Proposed Performance Index #2 


The second of the proposed performance metrics relies on the initial uncertainty 
in mean anomaly. It is given by; 




A 

= sup 






1/2 


(3.21) 


Using Equation (3.11) and (3.12), the theoretical behavior for Metric #2 can be 
expressed analytically in the following manner; 


ri"'{t^ = Inz sup 
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1/2 


(3.22) 


(3.23) 


(3.24) 


Given Equation (3.24) and points selected from the k<j ellipsoid, the maximum 
T]'” occurs when 5Mq = Q and |<7a*| is a maximum. Erom Equation (3.24), the 
theoretical rj"' can be expressed analytically as; 


V 


m 


= Itzz 


/ 

\ m 

^m+1 ( 



(3.25) 


The numerical results from Equation (3.21) are shown in Eigure 12 for 10 orbits. 
As with previous examples, the range of cr^^. is 1x10'"^ to 3x10'^ and is 0.01°. 
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l8t Order Ca$« 



2nd Order Case 



Figure 12. Proposed Performanee Metrie #2-10 Orbits 


Table 11 summarizes the results of the numerical versus theoretical comparison 
for Metric #2. As expected, the total change in the metric is larger for the first-order 
approximation. One can see there is good agreement between the numerical results 
obtained through Equation (3.21) and the predicted results from Equation (3.25). 
Interestingly, the behavior of Metric 2 is consistent throughout the range of uncertainty in 
So* and does not appear to be effected by the relative behavior of initial uncertainties in 
mean anomaly and semi-major axis. 



1.000x10'^ 

8.250x10'^ 

1.550x10'^ 

2.275x10'^ 

3.000x10'^ 

1 St Order Theoretical 

3.5343x10'^ 

0.2915 

0.5478 

0.8040 

1.0603 

1 St Order Numerical 

3.5355x10'^ 

0.2924 

0.5507 

0.8105 

1.0715 

2nd Order Theoretical 

1.2370x10'^ 

8.4193x10-'^ 

2.9719x10'^ 

6.4023x10'^ 

1.1133x10'^ 

2nd Order Numerical 

1.2374x10'^ 

8.4428x10-'^ 

2.9874x10'^ 

6.4518x10'^ 

1.1247x10'^ 


Table 11. Proposed Performance Metric #2-10 Orbits 
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4. Summary of Metric Performance 


The local nonlinearity index developed by Park and Scheeres was developed to 
decide the sufficient order of a higher order solution. The original metric as well as 
Index #1 represent the percent difference between the approximation and the true 
solution. This is very useful if one defines a threshold of accurateness for an 
approximation. Ultimately, the metrics do not give definitive information as to when the 
first-order approximation is no longer representing the system dynamics. 

D. RELATIVE BEHAVIOR OF INITIAL UNCERTAINTIES 

1. Overview 

Thus far it has been suggested that for the =0.0001 and =0.01° case, 
there is very good agreement over extended periods of propagation. This prompted 
further investigation into the relationship between the initial semi-major axis uncertainty 
and the mean anomaly. 

This section will examine the relative behavior between the initial uncertainties 
and the higher order effects on SM. Additionally, a performance metric is presented that 
predicts how long the linear approximation is valid given a target distribution and initial 
uncertainties. 

2. First-order Approximation at Small Initial Uncertainty 

Figure 13 shows the initial distribution of 10^ sample points propagated using the 
complete system dynamics, and then the 3cr uncertainty ellipse propagated using the first 
and second-order approximations (green ellipse and blue ellipse respectively). At 100 
orbits, there is very little difference between the first and second-order approximations 
and the first-order approximation still encapsulates nearly 99% of the initial distribution. 
The initial 3cr ellipse for the simulation captured 99.92% of the 10^ points. MATLAB’s 
pseudorandom number function randn() was used to generate the initial distribution. One 
can see that even after 100 orbits, the linear approximation identified by the green 3cr 
ellipse is representing the true system dynamics fairly well. The blue ellipse, which is 
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slightly distorted, is the seeond-order approximation, but overall there is very little 
difference. 


3ct Covariance After 100 Orbits 
Distribution Captured by 3a Ellipse = 0 9892 


'a 

>N 



• Monte Carlo - 1st Order - 2nd Order 


Figure 13. 3a Uncertainty Ellipse and Gaussian Distribution Propagated for 100 Orbits - 

= 0.7 km 


Likewise, Figure 14 shows the probability distribution curve using the same initial 
uncertainty in semi-major axis. At 100 orbits, there is near perfect agreement in hrst- 
order, second-order and exact propagations. Table 12 summarizes the results for the 
probability distribution out to 100 orbits. 
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Linear Propagation 2nd Order Propagation 




Exact Propagation 



k 


0 Orbits 
20 Orbits 
40 Orbits 
60 Orbits 
80 Orbits 
100 Orbits 
Theoretical 


Figure 14. Probability Distribution Comparison, 100 Orbits - <Jg^- 0.7 km 



Probability Distribution cr^^. = 0.0001 

Orbits 

2 

4 

6 

8 

10 

20 

40 

60 

80 

100 

Linear 

k @ 98.9% 

2.99 

2.99 

2.99 

3.00 

3.00 

3.00 

3.01 

3.00 

3.00 

3.02 

F(k) @ k=3> 

0.989 

0.989 

0.989 

0.989 

0.989 

0.989 

0.989 

0.989 

0.989 

0.988 

2nd 

Order 

k @ 98.9% 

2.99 

2.99 

2.99 

3.00 

3.00 

3.00 

3.01 

3.00 

3.00 

3.01 

F(k) @ k=3 

0.989 

0.989 

0.989 

0.989 

0.989 

0.989 

0.989 

0.989 

0.989 

0.988 

Exact 

k @ 98.9% 

3.02 

3.02 

3.02 

3.03 

3.03 

3.03 

3.03 

3.03 

3.03 

3.05 

F(k) @ k=3) 

0.988 

0.988 

0.988 

0.988 

0.988 

0.988 

0.988 

0.988 

0.987 

0.987 


Table 12. Probability Distribution - = 0.7 km 


Finally, looking at the results of the Horwood Prediction Error Cost Function, the 
Mahalanobis distance, and the performance indices #1 and #2, all of these metrics 
suggest that the cr^^. =0.0001 case, there is excellent agreement between first-order 
approximations. 

Initially, given that the equations of motion are a function of semi-major axis, it 


seems intuitive that with smaller initial uncertainty in semi-major axis, the first-order 
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approximation would remain valid for longer propagation time. However, the results at 
the cr^^. =0.0001 / (jgi^ =0.001° case show the first-order approximation captures the 
system dynamics out to extremely extended propagation times, past even 1000 orbits. 
Interestingly it was found that if the initial mean anomaly uncertainty was varied so that 
it was much smaller than the initial uncertainty in semi-major axis, the validity of the 
approximation was not as good at extended orbital periods. 

Figure 15 shows the probability distribution as well as the covariance propagation 
out to 100 orbits using cr^^. =0.0001. The uncertainty in mean anomaly however is 
changed from the original 0.01° to =0.001°. There is still fairly good agreement at 

100 orbits, but not as good as the =0.01° case. The right of Figure 15 shows 
approximately 92% of the distribution captured by the first-order approximation at 100 
orbits versus 98.9% at 100 orbits for =0.01°. This suggests that the validity of the 
first-order approximation is a function of the initial mean anomaly uncertainty, which 
was unexpected. 



Exact Propagation 
= 0 0001 
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80 Orbits 
100 Orbits 
Theoretical 


3a Covariance After 100 Orbits 
Distribution Captured by 3a Ellipse = 0.9182 
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Figure 15. Scr Uncertainty Ellipse and Gaussian Distribution Propagated for 100 Orbits - 

0.7 km and = 0.001° 
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3. 


First-order Approximation Error 


In this section, a comparison of the first-order approximation and the exact 
solution is made in order to obtain an expression for the validity of the linear 
approximation as a function of time and initial uncertainties. In doing so, this will also 
investigate the relative behavior between initial uncertainty in mean anomaly and semi¬ 
major axis and the influence on validity of propagation time. 

The exact solution for difference in mean anomaly expressed in non-dimensional 
terms is: 


= 2 ^ 


^ + 


-1 


t + SM, 


(3.26) 


The first-order approximation is: 


Order = 


(3.27) 


Define y as the first-order approximation error: 

T — Exact ~ ^^IstOrder 


(3.28) 


The difference, y is given by Equation (3.29): 


7 = 
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2/r 
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' ^ 3 1 

t + SMq 


[i 

(l + (^a*) 



- -3n5a T + 5 Mq (3.29) 


Simplified, Equation (3.29) becomes: 
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+ -Sa -1 
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l/TT 


(3.30) 
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It is notable that Equation (3.30) is not a function of This is inconsistent 

with previous results. In particular, Figures 13 through 15 suggest the initial uncertainty 
in 5Mq can have a significant influence on the first-order to exact solution error. 

Figure 16 shows y at 1000 orbits for a Sa range of ±60 km (60 km representing 
the approximate 3cr upper bound for Sa presented thus far in the thesis). A significant 
observation from Equation (3.30) and Figure 16 is that when 5a = Otheny = 0. The fact 
that the first-order approximation error is proportional to Sa* is a contributing factor to 
why the results thus far show good agreement for the Sa* = 0.0001. This is only part of 
the solution, as Equation (3.30) and Figure 16 fail to show any correlation between initial 
uncertainty in mean anomaly and the first-order approximation validity. 


0.9 
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6a* 

Figure 16. First-Order Approximation Error 

4, Linear Approximation Performance Metric 

As stated previously, the ultimate goal is to determine if there is a performance 
metric that can be used to identify when the nonlinearities in the equations of motion 


First-Order Approximation Error 

T I I I T 
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render the first-order approximation invalid. The following is an attempt to analytieally 
relate the validity of the first-order approximation with initial uneertainty in dM^. 

Given that previous results show the validity of the linear approximation is a 
funetion of the ka uneertainties were substituted into Equation (3.30) and y was set 
equal to k<jgf^^ in order to determine if this was the time when the linear approximation 
was no longer eapturing the system dynamies. 




+ — k(j . -1 
2 


27rT = [kcTg^^) 


(3.31) 


Equation (3.31) allows us to obtain a time in terms of number of orbits when the 
differenee between the first-order approximation and exaet solution is equal to 
This is given by Equation (3.32) 


r 








1 1 





In 


(3.32) 


Table 13 summarizes the results for A=3 using = 0.0T. The same range of 
initial uneertainties in semi-major axis used throughout this thesis is presented here. 
Additionally for referenee, the probability distribution F(k) at r = is given as well as 
the value of y at 1000 orbits. 



1 .000x10'^ 

8.250x10-^ 

1.550x10-^ 

2.275x10'^ 

3.000x10'^ 


494.0 

7.276 

2.067 

0.9617 

0.5545 

F(k) at r = Tj 

SMfi 

0.9681 

0.9681 

0.9681 

0.9681 

0.9681 

y at r=1000 orbits 

1.0599x10'^ 

7.1958x10'^ 

0.25336 

0.54443 

0.94434 


Table 13. Summary of Performanee of 
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As shown in Table 13, The first-order approximation F(k) is still at approximately 
96.8% of the distribution when is equal to y , which is still significant agreement. 

It follows that setting equal to y does not give the correct metric for assessing the 
validity of the linear approximation. However, previous results show the validity of the 
linear approximation is a function of so there still is value in examining Equation 


(3.32). 


To further illustrate the behavior of , Figure 17 shows Equation (3.32) 
plotted over a range of ±20 km in (7g^. The behavior of at (7g^ = 0 is infinite; 


however the y-axis is scaled to 1500 orbits in order to show the relevant behavior. If s 


clear from Equation (3.32) that with a larger numerator, the behavior of the curves 
approaching infinity will extend outward, essentially giving a greater number of orbits 
with larger values of cr^^. 
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Figure 17. with (7g^^ = 0.0T 


With respect to the results for Equation (3.32) listed in Table 13, there is 
interesting behavior with F(k). One can see that the number of orbits obtained change 
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with different cr^^, but the probability distribution, F(k) at r = does not ehange with 
different cr^^.. Using this property will allow us to determine an expression for the 
number of orbits as a funetion of <j^^. and that target a speeifie pereentage of the 
distribution we want the first-order approximation to eapture. 





^(l + Kv))’ 


+ - 




In 


(3.33) 


Given values for cr^^. and , the number of orbits for a target distribution ean 
be determined experimentally using either the error ellipse propagation or the probability 
distribution ealeulation. Equation (3.33) ean then be solved for the eoeffieient x. We ean 
then use Equation (3.33) to determine the approximate number of orbits the first-order 
approximation is good for with a given a target distribution, as a funetion of cr^^. and 
^SM„- ^ summary of the eoeffieient, x for various target distributions is listed in Table 
14. 


Target Distribution 

0.95 

0.9 

0.85 

0.8 

0.75 

0.7 

0.65 

Coefficient 

1.35020 

2.35425 

3.36235 

4.46964 

5.75911 

7.29757 

9.15385 


Table 14. Summary of Coeffieients used for Equation (3.33) 


Ultimately the relative performanee of the first-order approximation with respeet 
to is due to the pereentage differenee from the higher order effeets. The mean 
anomaly error, as well as any higher order eontributions to mean anomaly error, are 
solely a funetion of 5a . The first-order approximation for the uneertainty propagation 
has exeellent agreement with the seeond-order and exaet solution at = 0.0001 and 
agj^ = 0.0r beeause the higher order eontribution from 5a* is very small relative to (Jgj^j. 
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V. CONCLUSIONS 


A, KEY POINTS AND RECOMMENDATIONS 

In general, it was proven difficult to determine a definitive measure of when the 
covariance propagated using a linear approximation was no longer Gaussian. The 
cumulative distribution comparisons allow one to assess the performance of the 
approximation given initial uncertainty. The propagation of the covariance ellipsoid 
along with the Monte Carlo simulations provides a good method of visualizing the 
performance. Using these methods however, much like using the Mahalanobis distance, 
relies on the user to define a threshold that is no longer acceptable. 

The metrics such as the prediction error cost function and the local nonlinearity 
index provide a way of assessing the significant order of a higher order solution’s ability 
to fully capture a system’s dynamics. These metrics however did not provide any 
information on which initial uncertainties and their size cause propagated covariance to 
become non-Gaussian. 

The performance of the linear approximation with small initial uncertainty in cr^^. 
prompted investigation into the relative behavior between the initial uncertainties. This 
revealed that although the higher order effects are a function of Sa, they only have a 
significant impact when they are large in comparison to the uncertainty ellipse minor 
axis, which is a function of the initial uncertainty in mean anomaly. During the 
investigation into the relative behaviors, a linear approximation performance metric was 
developed. Although it is only applicable to the state as defined in this thesis, a similar 
approach could be used with a more complex model. 

B, AREAS OF FURTHER RESEARCH 

The system used in this thesis is very basic, consisting of only two-body motion 
in the semi-major axis/mean anomaly state space. Previous work referenced throughout 
this thesis has shown that different coordinate systems can produce vastly different 
results with respect to nonlinear behavior in first-order approximations. Further research 
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using different coordinate systems as well as using more complex models is warranted. 
Additionally, the linear approximation performance metric developed in this thesis could 
be modified to be made applicable to a real-world model and be used to predict more 
definitively how long a linear approximation is valid. 
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APPENDIX B. ADDITIONAL PLOTS 


Additional Cost Function Plots 



Figure 18. Cost Function Term 1 and 2 for 10 Orbits 
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Figure 19. Cost Function Term 1 and 2 for 100 Orbits 
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Figure 20. Total Cost Metric for 100 Orbits, First and Second-Order Approximation 

2, Additional Probability Distribution Plots 
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Figure 21. Probability Distribution Comparison, 10 Orbits - = 0.7 km 
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Figure 22. Probability Distribution Comparison, 10 Orbits - <Jg^ = 5.775 km 
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Figure 23. Probability Distribution Comparison, 10 Orbits - = 10.85 km 
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Figure 24. Probability Distribution Comparison, 100 Orbits - = 10.85 km 
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Figure 25. Probability Distribution Comparison, 10 Orbits - = 15.925 km 
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Figure 26. Probability Distribution Comparison, 100 Orbits - (Jg^ = 
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Figure 27. Probability Distribution Comparison, 10 Orbits - o-g 
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Figure 28. Probability Distribution Comparison, 100 Orbits - <Jg^ = 21 km 


3, Mean Anomaly Error and In-Track Error Plots 


1st Order Mean Anomaly Error for 10 Orbits 
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Figure 29. Mean Anomaly Error, 10 Orbits - Sa^ = 3(7^^., 5M^ - 0 
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1st Order Mean Anomaly Enor for 100 Orbits 
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Figure 30. Mean Anomaly Error, 100 Orbits - Sal - ~ ® 


1st Order In-Track Error for 10 Orbits 
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Figure 31. In-Traek Error, 10 Orbits - Sal ~ ® 
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1st Order In-Track Error for 100 Orbits 
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Figure 32. In-Track Error, 100 Orbits - Sal - ^ 
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